// Copyright (C) 2006, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#include <string>
#include <cassert>
#include <cfloat>
#include <cmath>
#include "CoinPragma.hpp"
#include "OsiSolverInterface.hpp"
#include "OsiAuxInfo.hpp"
#include "OsiSolverBranch.hpp"
#include "CoinWarmStartBasis.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinTime.hpp"
#include "CoinSort.hpp"
#include "CoinFinite.hpp"
#include "OsiChooseVariable.hpp"

OsiChooseVariable::OsiChooseVariable()
  : goodObjectiveValue_(COIN_DBL_MAX)
  , upChange_(0.0)
  , downChange_(0.0)
  , goodSolution_(NULL)
  , list_(NULL)
  , useful_(NULL)
  , solver_(NULL)
  , status_(-1)
  , bestObjectIndex_(-1)
  , bestWhichWay_(-1)
  , firstForcedObjectIndex_(-1)
  , firstForcedWhichWay_(-1)
  , numberUnsatisfied_(0)
  , numberStrong_(0)
  , numberOnList_(0)
  , numberStrongDone_(0)
  , numberStrongIterations_(0)
  , numberStrongFixed_(0)
  , trustStrongForBound_(true)
  , trustStrongForSolution_(true)
{
}

OsiChooseVariable::OsiChooseVariable(const OsiSolverInterface *solver)
  : goodObjectiveValue_(COIN_DBL_MAX)
  , upChange_(0.0)
  , downChange_(0.0)
  , goodSolution_(NULL)
  , solver_(solver)
  , status_(-1)
  , bestObjectIndex_(-1)
  , bestWhichWay_(-1)
  , firstForcedObjectIndex_(-1)
  , firstForcedWhichWay_(-1)
  , numberUnsatisfied_(0)
  , numberStrong_(0)
  , numberOnList_(0)
  , numberStrongDone_(0)
  , numberStrongIterations_(0)
  , numberStrongFixed_(0)
  , trustStrongForBound_(true)
  , trustStrongForSolution_(true)
{
  // create useful arrays
  int numberObjects = solver_->numberObjects();
  list_ = new int[numberObjects];
  useful_ = new double[numberObjects];
}

OsiChooseVariable::OsiChooseVariable(const OsiChooseVariable &rhs)
{
  goodObjectiveValue_ = rhs.goodObjectiveValue_;
  upChange_ = rhs.upChange_;
  downChange_ = rhs.downChange_;
  status_ = rhs.status_;
  bestObjectIndex_ = rhs.bestObjectIndex_;
  bestWhichWay_ = rhs.bestWhichWay_;
  firstForcedObjectIndex_ = rhs.firstForcedObjectIndex_;
  firstForcedWhichWay_ = rhs.firstForcedWhichWay_;
  numberUnsatisfied_ = rhs.numberUnsatisfied_;
  numberStrong_ = rhs.numberStrong_;
  numberOnList_ = rhs.numberOnList_;
  numberStrongDone_ = rhs.numberStrongDone_;
  numberStrongIterations_ = rhs.numberStrongIterations_;
  numberStrongFixed_ = rhs.numberStrongFixed_;
  trustStrongForBound_ = rhs.trustStrongForBound_;
  trustStrongForSolution_ = rhs.trustStrongForSolution_;
  solver_ = rhs.solver_;
  if (solver_) {
    int numberObjects = solver_->numberObjects();
    int numberColumns = solver_->getNumCols();
    if (rhs.goodSolution_) {
      goodSolution_ = CoinCopyOfArray(rhs.goodSolution_, numberColumns);
    } else {
      goodSolution_ = NULL;
    }
    list_ = CoinCopyOfArray(rhs.list_, numberObjects);
    useful_ = CoinCopyOfArray(rhs.useful_, numberObjects);
  } else {
    goodSolution_ = NULL;
    list_ = NULL;
    useful_ = NULL;
  }
}

OsiChooseVariable &
OsiChooseVariable::operator=(const OsiChooseVariable &rhs)
{
  if (this != &rhs) {
    delete[] goodSolution_;
    delete[] list_;
    delete[] useful_;
    goodObjectiveValue_ = rhs.goodObjectiveValue_;
    upChange_ = rhs.upChange_;
    downChange_ = rhs.downChange_;
    status_ = rhs.status_;
    bestObjectIndex_ = rhs.bestObjectIndex_;
    bestWhichWay_ = rhs.bestWhichWay_;
    firstForcedObjectIndex_ = rhs.firstForcedObjectIndex_;
    firstForcedWhichWay_ = rhs.firstForcedWhichWay_;
    numberUnsatisfied_ = rhs.numberUnsatisfied_;
    numberStrong_ = rhs.numberStrong_;
    numberOnList_ = rhs.numberOnList_;
    numberStrongDone_ = rhs.numberStrongDone_;
    numberStrongIterations_ = rhs.numberStrongIterations_;
    numberStrongFixed_ = rhs.numberStrongFixed_;
    trustStrongForBound_ = rhs.trustStrongForBound_;
    trustStrongForSolution_ = rhs.trustStrongForSolution_;
    solver_ = rhs.solver_;
    if (solver_) {
      int numberObjects = solver_->numberObjects();
      int numberColumns = solver_->getNumCols();
      if (rhs.goodSolution_) {
        goodSolution_ = CoinCopyOfArray(rhs.goodSolution_, numberColumns);
      } else {
        goodSolution_ = NULL;
      }
      list_ = CoinCopyOfArray(rhs.list_, numberObjects);
      useful_ = CoinCopyOfArray(rhs.useful_, numberObjects);
    } else {
      goodSolution_ = NULL;
      list_ = NULL;
      useful_ = NULL;
    }
  }
  return *this;
}

OsiChooseVariable::~OsiChooseVariable()
{
  delete[] goodSolution_;
  delete[] list_;
  delete[] useful_;
}

// Clone
OsiChooseVariable *
OsiChooseVariable::clone() const
{
  return new OsiChooseVariable(*this);
}
// Set solver and redo arrays
void OsiChooseVariable::setSolver(const OsiSolverInterface *solver)
{
  solver_ = solver;
  delete[] list_;
  delete[] useful_;
  // create useful arrays
  int numberObjects = solver_->numberObjects();
  list_ = new int[numberObjects];
  useful_ = new double[numberObjects];
}

// Initialize
int OsiChooseVariable::setupList(OsiBranchingInformation *info, bool initialize)
{
  if (initialize) {
    status_ = -2;
    delete[] goodSolution_;
    bestObjectIndex_ = -1;
    numberStrongDone_ = 0;
    numberStrongIterations_ = 0;
    numberStrongFixed_ = 0;
    goodSolution_ = NULL;
    goodObjectiveValue_ = COIN_DBL_MAX;
  }
  numberOnList_ = 0;
  numberUnsatisfied_ = 0;
  int numberObjects = solver_->numberObjects();
  assert(numberObjects);
  double check = 0.0;
  int checkIndex = 0;
  int bestPriority = COIN_INT_MAX;
  // pretend one strong even if none
  int maximumStrong = numberStrong_ ? std::min(numberStrong_, numberObjects) : 1;
  int putOther = numberObjects;
  int i;
  for (i = 0; i < maximumStrong; i++) {
    list_[i] = -1;
    useful_[i] = 0.0;
  }
  OsiObject **object = info->solver_->objects();
  // Say feasible
  bool feasible = true;
  for (i = 0; i < numberObjects; i++) {
    int way;
    double value = object[i]->infeasibility(info, way);
    if (value > 0.0) {
      numberUnsatisfied_++;
      if (value == COIN_DBL_MAX) {
        // infeasible
        feasible = false;
        break;
      }
      int priorityLevel = object[i]->priority();
      // Better priority? Flush choices.
      if (priorityLevel < bestPriority) {
        for (int j = 0; j < maximumStrong; j++) {
          if (list_[j] >= 0) {
            int iObject = list_[j];
            list_[j] = -1;
            useful_[j] = 0.0;
            list_[--putOther] = iObject;
          }
        }
        bestPriority = priorityLevel;
        check = 0.0;
      }
      if (priorityLevel == bestPriority) {
        if (value > check) {
          //add to list
          int iObject = list_[checkIndex];
          if (iObject >= 0)
            list_[--putOther] = iObject; // to end
          list_[checkIndex] = i;
          useful_[checkIndex] = value;
          // find worst
          check = COIN_DBL_MAX;
          for (int j = 0; j < maximumStrong; j++) {
            if (list_[j] >= 0) {
              if (useful_[j] < check) {
                check = useful_[j];
                checkIndex = j;
              }
            } else {
              check = 0.0;
              checkIndex = j;
              break;
            }
          }
        } else {
          // to end
          list_[--putOther] = i;
        }
      } else {
        // to end
        list_[--putOther] = i;
      }
    }
  }
  // Get list
  numberOnList_ = 0;
  if (feasible) {
    for (i = 0; i < maximumStrong; i++) {
      if (list_[i] >= 0) {
        list_[numberOnList_] = list_[i];
        useful_[numberOnList_++] = -useful_[i];
      }
    }
    if (numberOnList_) {
      // Sort
      CoinSort_2(useful_, useful_ + numberOnList_, list_);
      // move others
      i = numberOnList_;
      for (; putOther < numberObjects; putOther++)
        list_[i++] = list_[putOther];
      assert(i == numberUnsatisfied_);
      if (!numberStrong_)
        numberOnList_ = 0;
    }
  } else {
    // not feasible
    numberUnsatisfied_ = -1;
  }
  return numberUnsatisfied_;
}
/* Choose a variable
   Returns - 
   -1 Node is infeasible
   0  Normal termination - we have a candidate
   1  All looks satisfied - no candidate
   2  We can change the bound on a variable - but we also have a strong branching candidate
   3  We can change the bound on a variable - but we have a non-strong branching candidate
   4  We can change the bound on a variable - no other candidates
   We can pick up branch from whichObject() and whichWay()
   We can pick up a forced branch (can change bound) from whichForcedObject() and whichForcedWay()
   If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
*/
int OsiChooseVariable::chooseVariable(OsiSolverInterface *solver, OsiBranchingInformation *, bool)
{
  if (numberUnsatisfied_) {
    bestObjectIndex_ = list_[0];
    bestWhichWay_ = solver->object(bestObjectIndex_)->whichWay();
    firstForcedObjectIndex_ = -1;
    firstForcedWhichWay_ = -1;
    return 0;
  } else {
    return 1;
  }
}
// Returns true if solution looks feasible against given objects
bool OsiChooseVariable::feasibleSolution(const OsiBranchingInformation *info,
  const double *solution,
  int numberObjects,
  const OsiObject **objects)
{
  bool satisfied = true;
  const double *saveSolution = info->solution_;
  info->solution_ = solution;
  for (int i = 0; i < numberObjects; i++) {
    double value = objects[i]->checkInfeasibility(info);
    if (value > 0.0) {
      satisfied = false;
      break;
    }
  }
  info->solution_ = saveSolution;
  return satisfied;
}
// Saves a good solution
void OsiChooseVariable::saveSolution(const OsiSolverInterface *solver)
{
  delete[] goodSolution_;
  int numberColumns = solver->getNumCols();
  goodSolution_ = CoinCopyOfArray(solver->getColSolution(), numberColumns);
  goodObjectiveValue_ = solver->getObjSense() * solver->getObjValue();
}
// Clears out good solution after use
void OsiChooseVariable::clearGoodSolution()
{
  delete[] goodSolution_;
  goodSolution_ = NULL;
  goodObjectiveValue_ = COIN_DBL_MAX;
}

/*  This is a utility function which does strong branching on
    a list of objects and stores the results in OsiHotInfo.objects.
    On entry the object sequence is stored in the OsiHotInfo object
    and maybe more.
    It returns -
    -1 - one branch was infeasible both ways
     0 - all inspected - nothing can be fixed
     1 - all inspected - some can be fixed (returnCriterion==0)
     2 - may be returning early - one can be fixed (last one done) (returnCriterion==1) 
     3 - returning because max time
*/
int OsiChooseStrong::doStrongBranching(OsiSolverInterface *solver,
  OsiBranchingInformation *info,
  int numberToDo, int returnCriterion)
{

  // Might be faster to extend branch() to return bounds changed
  double *saveLower = NULL;
  double *saveUpper = NULL;
  int numberColumns = solver->getNumCols();
  solver->markHotStart();
  const double *lower = info->lower_;
  const double *upper = info->upper_;
  saveLower = CoinCopyOfArray(info->lower_, numberColumns);
  saveUpper = CoinCopyOfArray(info->upper_, numberColumns);
  numResults_ = 0;
  int returnCode = 0;
  double timeStart = CoinCpuTime();
  for (int iDo = 0; iDo < numberToDo; iDo++) {
    OsiHotInfo *result = results_ + iDo;
    // For now just 2 way
    OsiBranchingObject *branch = result->branchingObject();
    assert(branch->numberBranches() == 2);
    /*
      Try the first direction.  Each subsequent call to branch() performs the
      specified branch and advances the branch object state to the next branch
      alternative.)
    */
    OsiSolverInterface *thisSolver = solver;
    if (branch->boundBranch()) {
      // ordinary
      branch->branch(solver);
      // maybe we should check bounds for stupidities here?
      solver->solveFromHotStart();
    } else {
      // adding cuts or something
      thisSolver = solver->clone();
      branch->branch(thisSolver);
      // set hot start iterations
      int limit;
      thisSolver->getIntParam(OsiMaxNumIterationHotStart, limit);
      thisSolver->setIntParam(OsiMaxNumIteration, limit);
      thisSolver->resolve();
    }
    // can check if we got solution
    // status is 0 finished, 1 infeasible and 2 unfinished and 3 is solution
    int status0 = result->updateInformation(thisSolver, info, this);
    numberStrongIterations_ += thisSolver->getIterationCount();
    if (status0 == 3) {
      // new solution already saved
      if (trustStrongForSolution_) {
        info->cutoff_ = goodObjectiveValue_;
        status0 = 0;
      }
    }
    if (solver != thisSolver)
      delete thisSolver;
    // Restore bounds
    for (int j = 0; j < numberColumns; j++) {
      if (saveLower[j] != lower[j])
        solver->setColLower(j, saveLower[j]);
      if (saveUpper[j] != upper[j])
        solver->setColUpper(j, saveUpper[j]);
    }
    /*
      Try the next direction
    */
    thisSolver = solver;
    if (branch->boundBranch()) {
      // ordinary
      branch->branch(solver);
      // maybe we should check bounds for stupidities here?
      solver->solveFromHotStart();
    } else {
      // adding cuts or something
      thisSolver = solver->clone();
      branch->branch(thisSolver);
      // set hot start iterations
      int limit;
      thisSolver->getIntParam(OsiMaxNumIterationHotStart, limit);
      thisSolver->setIntParam(OsiMaxNumIteration, limit);
      thisSolver->resolve();
    }
    // can check if we got solution
    // status is 0 finished, 1 infeasible and 2 unfinished and 3 is solution
    int status1 = result->updateInformation(thisSolver, info, this);
    numberStrongDone_++;
    numberStrongIterations_ += thisSolver->getIterationCount();
    if (status1 == 3) {
      // new solution already saved
      if (trustStrongForSolution_) {
        info->cutoff_ = goodObjectiveValue_;
        status1 = 0;
      }
    }
    if (solver != thisSolver)
      delete thisSolver;
    // Restore bounds
    for (int j = 0; j < numberColumns; j++) {
      if (saveLower[j] != lower[j])
        solver->setColLower(j, saveLower[j]);
      if (saveUpper[j] != upper[j])
        solver->setColUpper(j, saveUpper[j]);
    }
    /*
      End of evaluation for this candidate variable. Possibilities are:
      * Both sides below cutoff; this variable is a candidate for branching.
      * Both sides infeasible or above the objective cutoff: no further action
      here. Break from the evaluation loop and assume the node will be purged
      by the caller.
      * One side below cutoff: Install the branch (i.e., fix the variable). Possibly break
      from the evaluation loop and assume the node will be reoptimised by the
      caller.
    */
    numResults_++;
    if (status0 == 1 && status1 == 1) {
      // infeasible
      returnCode = -1;
      break; // exit loop
    } else if (status0 == 1 || status1 == 1) {
      numberStrongFixed_++;
      if (!returnCriterion) {
        returnCode = 1;
      } else {
        returnCode = 2;
        break;
      }
    }
    bool hitMaxTime = (CoinCpuTime() - timeStart > info->timeRemaining_);
    if (hitMaxTime) {
      returnCode = 3;
      break;
    }
  }
  delete[] saveLower;
  delete[] saveUpper;
  // Delete the snapshot
  solver->unmarkHotStart();
  return returnCode;
}

// Given a candidate fill in useful information e.g. estimates
void OsiChooseVariable::updateInformation(const OsiBranchingInformation *info,
  int, OsiHotInfo *hotInfo)
{
  int index = hotInfo->whichObject();
  assert(index < solver_->numberObjects());
  //assert (branch<2);
  OsiObject **object = info->solver_->objects();
  upChange_ = object[index]->upEstimate();
  downChange_ = object[index]->downEstimate();
}
#if 1
// Given a branch fill in useful information e.g. estimates
void OsiChooseVariable::updateInformation(int index, int branch,
  double, double,
  int)
{
  assert(index < solver_->numberObjects());
  assert(branch < 2);
  OsiObject **object = solver_->objects();
  if (branch)
    upChange_ = object[index]->upEstimate();
  else
    downChange_ = object[index]->downEstimate();
}
#endif

//##############################################################################

void OsiPseudoCosts::gutsOfDelete()
{
  if (numberObjects_ > 0) {
    numberObjects_ = 0;
    numberBeforeTrusted_ = 0;
    delete[] upTotalChange_;
    upTotalChange_ = NULL;
    delete[] downTotalChange_;
    downTotalChange_ = NULL;
    delete[] upNumber_;
    upNumber_ = NULL;
    delete[] downNumber_;
    downNumber_ = NULL;
  }
}

void OsiPseudoCosts::gutsOfCopy(const OsiPseudoCosts &rhs)
{
  numberObjects_ = rhs.numberObjects_;
  numberBeforeTrusted_ = rhs.numberBeforeTrusted_;
  if (numberObjects_ > 0) {
    upTotalChange_ = CoinCopyOfArray(rhs.upTotalChange_, numberObjects_);
    downTotalChange_ = CoinCopyOfArray(rhs.downTotalChange_, numberObjects_);
    upNumber_ = CoinCopyOfArray(rhs.upNumber_, numberObjects_);
    downNumber_ = CoinCopyOfArray(rhs.downNumber_, numberObjects_);
  }
}

OsiPseudoCosts::OsiPseudoCosts()
  : upTotalChange_(NULL)
  , downTotalChange_(NULL)
  , upNumber_(NULL)
  , downNumber_(NULL)
  , numberObjects_(0)
  , numberBeforeTrusted_(0)
{
}

OsiPseudoCosts::~OsiPseudoCosts()
{
  gutsOfDelete();
}

OsiPseudoCosts::OsiPseudoCosts(const OsiPseudoCosts &rhs)
  : upTotalChange_(NULL)
  , downTotalChange_(NULL)
  , upNumber_(NULL)
  , downNumber_(NULL)
  , numberObjects_(0)
  , numberBeforeTrusted_(0)
{
  gutsOfCopy(rhs);
}

OsiPseudoCosts &
OsiPseudoCosts::operator=(const OsiPseudoCosts &rhs)
{
  if (this != &rhs) {
    gutsOfDelete();
    gutsOfCopy(rhs);
  }
  return *this;
}

void OsiPseudoCosts::initialize(int n)
{
  gutsOfDelete();
  numberObjects_ = n;
  if (numberObjects_ > 0) {
    upTotalChange_ = new double[numberObjects_];
    downTotalChange_ = new double[numberObjects_];
    upNumber_ = new int[numberObjects_];
    downNumber_ = new int[numberObjects_];
    CoinZeroN(upTotalChange_, numberObjects_);
    CoinZeroN(downTotalChange_, numberObjects_);
    CoinZeroN(upNumber_, numberObjects_);
    CoinZeroN(downNumber_, numberObjects_);
  }
}

//##############################################################################

OsiChooseStrong::OsiChooseStrong()
  : OsiChooseVariable()
  , shadowPriceMode_(0)
  , pseudoCosts_()
  , results_(NULL)
  , numResults_(0)
{
}

OsiChooseStrong::OsiChooseStrong(const OsiSolverInterface *solver)
  : OsiChooseVariable(solver)
  , shadowPriceMode_(0)
  , pseudoCosts_()
  , results_(NULL)
  , numResults_(0)
{
  // create useful arrays
  pseudoCosts_.initialize(solver_->numberObjects());
}

OsiChooseStrong::OsiChooseStrong(const OsiChooseStrong &rhs)
  : OsiChooseVariable(rhs)
  , shadowPriceMode_(rhs.shadowPriceMode_)
  , pseudoCosts_(rhs.pseudoCosts_)
  , results_(NULL)
  , numResults_(0)
{
}

OsiChooseStrong &
OsiChooseStrong::operator=(const OsiChooseStrong &rhs)
{
  if (this != &rhs) {
    OsiChooseVariable::operator=(rhs);
    shadowPriceMode_ = rhs.shadowPriceMode_;
    pseudoCosts_ = rhs.pseudoCosts_;
    delete[] results_;
    results_ = NULL;
    numResults_ = 0;
  }
  return *this;
}

OsiChooseStrong::~OsiChooseStrong()
{
  delete[] results_;
}

// Clone
OsiChooseVariable *
OsiChooseStrong::clone() const
{
  return new OsiChooseStrong(*this);
}
#define MAXMIN_CRITERION 0.85
// Initialize
int OsiChooseStrong::setupList(OsiBranchingInformation *info, bool initialize)
{
  if (initialize) {
    status_ = -2;
    delete[] goodSolution_;
    bestObjectIndex_ = -1;
    numberStrongDone_ = 0;
    numberStrongIterations_ = 0;
    numberStrongFixed_ = 0;
    goodSolution_ = NULL;
    goodObjectiveValue_ = COIN_DBL_MAX;
  }
  numberOnList_ = 0;
  numberUnsatisfied_ = 0;
  int numberObjects = solver_->numberObjects();
  if (numberObjects > pseudoCosts_.numberObjects()) {
    // redo useful arrays
    pseudoCosts_.initialize(numberObjects);
  }
  double check = -COIN_DBL_MAX;
  int checkIndex = 0;
  int bestPriority = COIN_INT_MAX;
  int maximumStrong = std::min(numberStrong_, numberObjects);
  int putOther = numberObjects;
  int i;
  for (i = 0; i < numberObjects; i++) {
    list_[i] = -1;
    useful_[i] = 0.0;
  }
  OsiObject **object = info->solver_->objects();
  // Get average pseudo costs and see if pseudo shadow prices possible
  int shadowPossible = shadowPriceMode_;
  if (shadowPossible) {
    for (i = 0; i < numberObjects; i++) {
      if (!object[i]->canHandleShadowPrices()) {
        shadowPossible = 0;
        break;
      }
    }
    if (shadowPossible) {
      int numberRows = solver_->getNumRows();
      const double *pi = info->pi_;
      double sumPi = 0.0;
      for (i = 0; i < numberRows; i++)
        sumPi += fabs(pi[i]);
      sumPi /= static_cast< double >(numberRows);
      // and scale back
      sumPi *= 0.01;
      info->defaultDual_ = sumPi; // switch on
      int numberColumns = solver_->getNumCols();
      int size = std::max(numberColumns, 2 * numberRows);
      info->usefulRegion_ = new double[size];
      CoinZeroN(info->usefulRegion_, size);
      info->indexRegion_ = new int[size];
    }
  }
  double sumUp = 0.0;
  double numberUp = 0.0;
  double sumDown = 0.0;
  double numberDown = 0.0;
  const double *upTotalChange = pseudoCosts_.upTotalChange();
  const double *downTotalChange = pseudoCosts_.downTotalChange();
  const int *upNumber = pseudoCosts_.upNumber();
  const int *downNumber = pseudoCosts_.downNumber();
  const int numberBeforeTrusted = pseudoCosts_.numberBeforeTrusted();
  for (i = 0; i < numberObjects; i++) {
    sumUp += upTotalChange[i];
    numberUp += upNumber[i];
    sumDown += downTotalChange[i];
    numberDown += downNumber[i];
  }
  double upMultiplier = (1.0 + sumUp) / (1.0 + numberUp);
  double downMultiplier = (1.0 + sumDown) / (1.0 + numberDown);
  // Say feasible
  bool feasible = true;
#if 0
  int pri[]={10,1000,10000};
  int priCount[]={0,0,0};
#endif
  for (i = 0; i < numberObjects; i++) {
    int way;
    double value = object[i]->infeasibility(info, way);
    if (value > 0.0) {
      numberUnsatisfied_++;
      if (value == COIN_DBL_MAX) {
        // infeasible
        feasible = false;
        break;
      }
      int priorityLevel = object[i]->priority();
#if 0
      for (int k=0;k<3;k++) {
	if (priorityLevel==pri[k])
	  priCount[k]++;
      }
#endif
      // Better priority? Flush choices.
      if (priorityLevel < bestPriority) {
        for (int j = maximumStrong - 1; j >= 0; j--) {
          if (list_[j] >= 0) {
            int iObject = list_[j];
            list_[j] = -1;
            useful_[j] = 0.0;
            list_[--putOther] = iObject;
          }
        }
        maximumStrong = std::min(maximumStrong, putOther);
        bestPriority = priorityLevel;
        check = -COIN_DBL_MAX;
        checkIndex = 0;
      }
      if (priorityLevel == bestPriority) {
        // Modify value
        sumUp = upTotalChange[i] + 1.0e-30;
        numberUp = upNumber[i];
        sumDown = downTotalChange[i] + 1.0e-30;
        numberDown = downNumber[i];
        double upEstimate = object[i]->upEstimate();
        double downEstimate = object[i]->downEstimate();
        if (shadowPossible < 2) {
          upEstimate = numberUp ? ((upEstimate * sumUp) / numberUp) : (upEstimate * upMultiplier);
          if (numberUp < numberBeforeTrusted)
            upEstimate *= (numberBeforeTrusted + 1.0) / (numberUp + 1.0);
          downEstimate = numberDown ? ((downEstimate * sumDown) / numberDown) : (downEstimate * downMultiplier);
          if (numberDown < numberBeforeTrusted)
            downEstimate *= (numberBeforeTrusted + 1.0) / (numberDown + 1.0);
        } else {
          // use shadow prices always
        }
        value = MAXMIN_CRITERION * std::min(upEstimate, downEstimate) + (1.0 - MAXMIN_CRITERION) * std::max(upEstimate, downEstimate);
        if (value > check) {
          //add to list
          int iObject = list_[checkIndex];
          if (iObject >= 0) {
            assert(list_[putOther - 1] < 0);
            list_[--putOther] = iObject; // to end
          }
          list_[checkIndex] = i;
          assert(checkIndex < putOther);
          useful_[checkIndex] = value;
          // find worst
          check = COIN_DBL_MAX;
          maximumStrong = std::min(maximumStrong, putOther);
          for (int j = 0; j < maximumStrong; j++) {
            if (list_[j] >= 0) {
              if (useful_[j] < check) {
                check = useful_[j];
                checkIndex = j;
              }
            } else {
              check = 0.0;
              checkIndex = j;
              break;
            }
          }
        } else {
          // to end
          assert(list_[putOther - 1] < 0);
          list_[--putOther] = i;
          maximumStrong = std::min(maximumStrong, putOther);
        }
      } else {
        // worse priority
        // to end
        assert(list_[putOther - 1] < 0);
        list_[--putOther] = i;
        maximumStrong = std::min(maximumStrong, putOther);
      }
    }
  }
#if 0
  printf("%d at %d, %d at %d and %d at %d\n",priCount[0],pri[0],
	 priCount[1],pri[1],priCount[2],pri[2]);
#endif
  // Get list
  numberOnList_ = 0;
  if (feasible) {
    for (i = 0; i < std::min(maximumStrong, putOther); i++) {
      if (list_[i] >= 0) {
        list_[numberOnList_] = list_[i];
        useful_[numberOnList_++] = -useful_[i];
      }
    }
    if (numberOnList_) {
      // Sort
      CoinSort_2(useful_, useful_ + numberOnList_, list_);
      // move others
      i = numberOnList_;
      for (; putOther < numberObjects; putOther++)
        list_[i++] = list_[putOther];
      assert(i == numberUnsatisfied_);
      if (!numberStrong_)
        numberOnList_ = 0;
    }
  } else {
    // not feasible
    numberUnsatisfied_ = -1;
  }
  // Get rid of any shadow prices info
  info->defaultDual_ = -1.0; // switch off
  delete[] info->usefulRegion_;
  delete[] info->indexRegion_;
  return numberUnsatisfied_;
}

void OsiChooseStrong::resetResults(int num)
{
  delete[] results_;
  numResults_ = 0;
  results_ = new OsiHotInfo[num];
}

/* Choose a variable
   Returns - 
   -1 Node is infeasible
   0  Normal termination - we have a candidate
   1  All looks satisfied - no candidate
   2  We can change the bound on a variable - but we also have a strong branching candidate
   3  We can change the bound on a variable - but we have a non-strong branching candidate
   4  We can change the bound on a variable - no other candidates
   We can pick up branch from whichObject() and whichWay()
   We can pick up a forced branch (can change bound) from whichForcedObject() and whichForcedWay()
   If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
*/
int OsiChooseStrong::chooseVariable(OsiSolverInterface *solver, OsiBranchingInformation *info, bool fixVariables)
{
  if (numberUnsatisfied_) {
    const double *upTotalChange = pseudoCosts_.upTotalChange();
    const double *downTotalChange = pseudoCosts_.downTotalChange();
    const int *upNumber = pseudoCosts_.upNumber();
    const int *downNumber = pseudoCosts_.downNumber();
    int numberBeforeTrusted = pseudoCosts_.numberBeforeTrusted();
    // Somehow we can get here with it 0 !
    if (!numberBeforeTrusted) {
      numberBeforeTrusted = 5;
      pseudoCosts_.setNumberBeforeTrusted(numberBeforeTrusted);
    }

    int numberLeft = std::min(numberStrong_ - numberStrongDone_, numberUnsatisfied_);
    int numberToDo = 0;
    resetResults(numberLeft);
    int returnCode = 0;
    bestObjectIndex_ = -1;
    bestWhichWay_ = -1;
    firstForcedObjectIndex_ = -1;
    firstForcedWhichWay_ = -1;
    double bestTrusted = -COIN_DBL_MAX;
    for (int i = 0; i < numberLeft; i++) {
      int iObject = list_[i];
      if (upNumber[iObject] < numberBeforeTrusted || downNumber[iObject] < numberBeforeTrusted) {
        results_[numberToDo++] = OsiHotInfo(solver, info,
          solver->objects(), iObject);
      } else {
        const OsiObject *obj = solver->object(iObject);
        double upEstimate = (upTotalChange[iObject] * obj->upEstimate()) / upNumber[iObject];
        double downEstimate = (downTotalChange[iObject] * obj->downEstimate()) / downNumber[iObject];
        double value = MAXMIN_CRITERION * std::min(upEstimate, downEstimate) + (1.0 - MAXMIN_CRITERION) * std::max(upEstimate, downEstimate);
        if (value > bestTrusted) {
          bestObjectIndex_ = iObject;
          bestWhichWay_ = upEstimate > downEstimate ? 0 : 1;
          bestTrusted = value;
        }
      }
    }
    int numberFixed = 0;
    if (numberToDo) {
      returnCode = doStrongBranching(solver, info, numberToDo, 1);
      if (returnCode >= 0 && returnCode <= 2) {
        if (returnCode) {
          returnCode = 4;
          if (bestObjectIndex_ >= 0)
            returnCode = 3;
        }
        for (int i = 0; i < numResults_; i++) {
          int iObject = results_[i].whichObject();
          double upEstimate;
          if (results_[i].upStatus() != 1) {
            assert(results_[i].upStatus() >= 0);
            upEstimate = results_[i].upChange();
          } else {
            // infeasible - just say expensive
            if (info->cutoff_ < 1.0e50)
              upEstimate = 2.0 * (info->cutoff_ - info->objectiveValue_);
            else
              upEstimate = 2.0 * fabs(info->objectiveValue_);
            if (firstForcedObjectIndex_ < 0) {
              firstForcedObjectIndex_ = iObject;
              firstForcedWhichWay_ = 0;
            }
            numberFixed++;
            if (fixVariables) {
              const OsiObject *obj = solver->object(iObject);
              OsiBranchingObject *branch = obj->createBranch(solver, info, 0);
              branch->branch(solver);
              delete branch;
            }
          }
          double downEstimate;
          if (results_[i].downStatus() != 1) {
            assert(results_[i].downStatus() >= 0);
            downEstimate = results_[i].downChange();
          } else {
            // infeasible - just say expensive
            if (info->cutoff_ < 1.0e50)
              downEstimate = 2.0 * (info->cutoff_ - info->objectiveValue_);
            else
              downEstimate = 2.0 * fabs(info->objectiveValue_);
            if (firstForcedObjectIndex_ < 0) {
              firstForcedObjectIndex_ = iObject;
              firstForcedWhichWay_ = 1;
            }
            numberFixed++;
            if (fixVariables) {
              const OsiObject *obj = solver->object(iObject);
              OsiBranchingObject *branch = obj->createBranch(solver, info, 1);
              branch->branch(solver);
              delete branch;
            }
          }
          double value = MAXMIN_CRITERION * std::min(upEstimate, downEstimate) + (1.0 - MAXMIN_CRITERION) * std::max(upEstimate, downEstimate);
          if (value > bestTrusted) {
            bestTrusted = value;
            bestObjectIndex_ = iObject;
            bestWhichWay_ = upEstimate > downEstimate ? 0 : 1;
            // but override if there is a preferred way
            const OsiObject *obj = solver->object(iObject);
            if (obj->preferredWay() >= 0 && obj->infeasibility())
              bestWhichWay_ = obj->preferredWay();
            if (returnCode)
              returnCode = 2;
          }
        }
      } else if (returnCode == 3) {
        // max time - just choose one
        bestObjectIndex_ = list_[0];
        bestWhichWay_ = 0;
        returnCode = 0;
      }
    } else {
      bestObjectIndex_ = list_[0];
    }
    if (bestObjectIndex_ >= 0) {
      OsiObject *obj = solver->objects()[bestObjectIndex_];
      obj->setWhichWay(bestWhichWay_);
    }
    if (numberFixed == numberUnsatisfied_ && numberFixed)
      returnCode = 4;
    return returnCode;
  } else {
    return 1;
  }
}
// Given a candidate  fill in useful information e.g. estimates
void OsiPseudoCosts::updateInformation(const OsiBranchingInformation *info,
  int branch, OsiHotInfo *hotInfo)
{
  int index = hotInfo->whichObject();
  assert(index < info->solver_->numberObjects());
  const OsiObject *object = info->solver_->object(index);
  assert(object->upEstimate() > 0.0 && object->downEstimate() > 0.0);
  assert(branch < 2);
  if (branch) {
    if (hotInfo->upStatus() != 1) {
      assert(hotInfo->upStatus() >= 0);
      upTotalChange_[index] += hotInfo->upChange() / object->upEstimate();
      upNumber_[index]++;
    } else {
#if 0
      // infeasible - just say expensive
      if (info->cutoff_<1.0e50)
	upTotalChange_[index] += 2.0*(info->cutoff_-info->objectiveValue_)/object->upEstimate();
      else
	upTotalChange_[index] += 2.0*fabs(info->objectiveValue_)/object->upEstimate();
#endif
    }
  } else {
    if (hotInfo->downStatus() != 1) {
      assert(hotInfo->downStatus() >= 0);
      downTotalChange_[index] += hotInfo->downChange() / object->downEstimate();
      downNumber_[index]++;
    } else {
#if 0
      // infeasible - just say expensive
      if (info->cutoff_<1.0e50)
	downTotalChange_[index] += 2.0*(info->cutoff_-info->objectiveValue_)/object->downEstimate();
      else
	downTotalChange_[index] += 2.0*fabs(info->objectiveValue_)/object->downEstimate();
#endif
    }
  }
}
#if 1
// Given a branch fill in useful information e.g. estimates
void OsiPseudoCosts::updateInformation(int index, int branch,
  double changeInObjective,
  double changeInValue,
  int status)
{
  //assert (index<solver_->numberObjects());
  assert(branch < 2);
  assert(changeInValue > 0.0);
  assert(branch < 2);
  if (branch) {
    if (status != 1) {
      assert(status >= 0);
      upTotalChange_[index] += changeInObjective / changeInValue;
      upNumber_[index]++;
    }
  } else {
    if (status != 1) {
      assert(status >= 0);
      downTotalChange_[index] += changeInObjective / changeInValue;
      downNumber_[index]++;
    }
  }
}
#endif

OsiHotInfo::OsiHotInfo()
  : originalObjectiveValue_(COIN_DBL_MAX)
  , changes_(NULL)
  , iterationCounts_(NULL)
  , statuses_(NULL)
  , branchingObject_(NULL)
  , whichObject_(-1)
{
}

OsiHotInfo::OsiHotInfo(OsiSolverInterface *solver,
  const OsiBranchingInformation *info,
  const OsiObject *const *objects,
  int whichObject)
  : originalObjectiveValue_(COIN_DBL_MAX)
  , whichObject_(whichObject)
{
  originalObjectiveValue_ = info->objectiveValue_;
  const OsiObject *object = objects[whichObject_];
  // create object - "down" first
  branchingObject_ = object->createBranch(solver, info, 0);
  // create arrays
  int numberBranches = branchingObject_->numberBranches();
  changes_ = new double[numberBranches];
  iterationCounts_ = new int[numberBranches];
  statuses_ = new int[numberBranches];
  CoinZeroN(changes_, numberBranches);
  CoinZeroN(iterationCounts_, numberBranches);
  CoinFillN(statuses_, numberBranches, -1);
}

OsiHotInfo::OsiHotInfo(const OsiHotInfo &rhs)
{
  originalObjectiveValue_ = rhs.originalObjectiveValue_;
  whichObject_ = rhs.whichObject_;
  if (rhs.branchingObject_) {
    branchingObject_ = rhs.branchingObject_->clone();
    int numberBranches = branchingObject_->numberBranches();
    changes_ = CoinCopyOfArray(rhs.changes_, numberBranches);
    iterationCounts_ = CoinCopyOfArray(rhs.iterationCounts_, numberBranches);
    statuses_ = CoinCopyOfArray(rhs.statuses_, numberBranches);
  } else {
    branchingObject_ = NULL;
    changes_ = NULL;
    iterationCounts_ = NULL;
    statuses_ = NULL;
  }
}

OsiHotInfo &
OsiHotInfo::operator=(const OsiHotInfo &rhs)
{
  if (this != &rhs) {
    delete branchingObject_;
    delete[] changes_;
    delete[] iterationCounts_;
    delete[] statuses_;
    originalObjectiveValue_ = rhs.originalObjectiveValue_;
    whichObject_ = rhs.whichObject_;
    if (rhs.branchingObject_) {
      branchingObject_ = rhs.branchingObject_->clone();
      int numberBranches = branchingObject_->numberBranches();
      changes_ = CoinCopyOfArray(rhs.changes_, numberBranches);
      iterationCounts_ = CoinCopyOfArray(rhs.iterationCounts_, numberBranches);
      statuses_ = CoinCopyOfArray(rhs.statuses_, numberBranches);
    } else {
      branchingObject_ = NULL;
      changes_ = NULL;
      iterationCounts_ = NULL;
      statuses_ = NULL;
    }
  }
  return *this;
}

OsiHotInfo::~OsiHotInfo()
{
  delete branchingObject_;
  delete[] changes_;
  delete[] iterationCounts_;
  delete[] statuses_;
}

// Clone
OsiHotInfo *
OsiHotInfo::clone() const
{
  return new OsiHotInfo(*this);
}
/* Fill in useful information after strong branch 
 */
int OsiHotInfo::updateInformation(const OsiSolverInterface *solver, const OsiBranchingInformation *info,
  OsiChooseVariable *choose)
{
  int iBranch = branchingObject_->branchIndex() - 1;
  assert(iBranch >= 0 && iBranch < branchingObject_->numberBranches());
  iterationCounts_[iBranch] += solver->getIterationCount();
  int status;
  if (solver->isProvenOptimal())
    status = 0; // optimal
  else if (solver->isIterationLimitReached()
    && !solver->isDualObjectiveLimitReached())
    status = 2; // unknown
  else
    status = 1; // infeasible
  // Could do something different if we can't trust
  double newObjectiveValue = solver->getObjSense() * solver->getObjValue();
  changes_[iBranch] = std::max(0.0, newObjectiveValue - originalObjectiveValue_);
  // we might have got here by primal
  if (choose->trustStrongForBound()) {
    if (!status && newObjectiveValue >= info->cutoff_) {
      status = 1; // infeasible
      changes_[iBranch] = 1.0e100;
    }
  }
  statuses_[iBranch] = status;
  if (!status && choose->trustStrongForSolution() && newObjectiveValue < choose->goodObjectiveValue()) {
    // check if solution
    const OsiSolverInterface *saveSolver = info->solver_;
    info->solver_ = solver;
    const double *saveLower = info->lower_;
    info->lower_ = solver->getColLower();
    const double *saveUpper = info->upper_;
    info->upper_ = solver->getColUpper();
    // also need to make sure bounds OK as may not be info solver
#if 0
    if (saveSolver->getMatrixByCol()) {
	const CoinBigIndex * columnStart = info->columnStart_;
	assert (saveSolver->getMatrixByCol()->getVectorStarts()==columnStart);
    }
#endif
    if (choose->feasibleSolution(info, solver->getColSolution(), solver->numberObjects(),
          const_cast< const OsiObject ** >(solver->objects()))) {
      // put solution somewhere
      choose->saveSolution(solver);
      status = 3;
    }
    info->solver_ = saveSolver;
    info->lower_ = saveLower;
    info->upper_ = saveUpper;
  }
  // Now update - possible strong branching info
  choose->updateInformation(info, iBranch, this);
  return status;
}

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/
